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FOREWORD 


This  research  report  is  concerned  with  the  dynamic  response  of  unidirectional 
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ture  and  Solid  Mechanics  at  Lehigh  University.  The  Principal  Investigator  of 
the  project  is  Professor  George  C.  Sih  and  the  Associate  Investigator  is  Dr.  E. 

P.  Chen  who  has  since  left  Lehigh  University  and  joined  the  Sandia  Laboratory 
in  New  Mexico.  The  authors  are  grateful  to  the  NASA  Project  Manager,  Dr.  Christos 
C.  Chamis  who  has  carefully  reviewed  this  report  and  provided  a  number  of  concrete 
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OFF-AXIS  IMPACT  OF  UNIDIRECTIONAL  COMPOSITES  WITH  CRACKS: 
DYNAMIC  STRESS  INTENSIFICATION 


by 

G.  C.  Sih 

Institute  of  Fracture  and  Solid  Mechanics 
Lehigh  University 
Bethlehem,  Pennsylvania  18015 

and 

E.  P.  Chen 
Sandia  Laboratories 
Albuquerque,  New  Mexico  87115 

ABSTRACT 

The  dynamic  response  of  unidirectional  composites  under  off-axis  (angle 
loading)  impact  is  analyzed  by  assuming  that  the  composite  contains  an  initial 
flaw  in  the  matrix  material.  Because  of  the  complexities  that  arise  from  the 
interaction  of  waves  scattered  by  the  crack  with  those  reflected  by  the  inter¬ 
faces  within  the  composite,  dynamic  analyses  of  composites  with  cracks  have  been 
treated  only  for  a  few  simple  cases.  One  of  the  objectives  of  the  present  work 
is  to  develop  an  effective  analytical  method  for  determining  dynamic  stress  so¬ 
lutions.  This  will  not  only  lead  to  an  in-depth  understanding  of  the  failure 
of  composites  due  to  impact  but  also  provide  reliable  solutions  that  can  guide 
the  development  of  numerical  methods. 

The  analysis  method  utilizes  Fourier  transform  for  the  space  variable  and 
Laplace  transform  for  the  time  variable.  The  time-dependent  angle  loading  is 


*This  work  was  completed  during  Dr.  Chen's  tenure  at  Lehigh  University. 
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separated  into  two  parts:  one  being  symmetric  and  the  other  skew- symmetric  with 
reference  to  the  crack  plane.  By  means  of  superposition,  the  transient  boundary 
conditions  consist  of  applying  normal  and  shear  tractions  to  a  crack  embedded  in 
the  matrix  of  the  unidirectional  composite.  Mathematically,  these  conditions 
reduce  the  problem  to  a  system  of  dual  integral  equations  which  are  solved  in 
the  Laplace  transform  plane  for  the  transform  of  the  dynamic  stress  intensity 
factor.  The  time  inversion  is  carried  out  numerically  for  various  combinations 
of  the  material  properties  of  the  composite  and  the  results  are  displayed  graph¬ 
ically. 


INTRODUCTION 

Past  work  on  the  development  of  high  performance  composite  materials  was 
mainly  concerned  with  achieving  high  strength  and  modulus.  This  requirement 
alone,  however,  may  result  in  a  composite  that  is  excessively  brittle  and  lacks 
the  ability  to  resist  impact  loading.  The  energy  absorption  or  toughness  of  the 
composite  is  also  an  important  property  that  must  be  accounted  for  in  addition 
to  strength  and  stiffness. 

The  concept  of  fracture  toughness  has  mostly  been  applied  to  homogeneous 
isotropic  materials  [1]  based  on  the  linear  fracture  mechanics  theories  such  as 
those  advanced  by  Griffith,  Irwin  and  others.  These  theories,  developed  for 
single-phase  materials,  have  had  limited  success  in  characterizing  the  fracture 
behavior  of  composites  which  are  inherently  nonhomogeneous  and  anisotropic. 

This  is  mainly  because  the  fracture  modes  in  composites  are  multi-facet  and  can 
include  interface  failure,  fiber  breaking,  matrix  fracture,  etc.  The  individ¬ 
ual  contribution  of  each  of  these  failure  modes  is  not  clearly  accounted  for 
and/or  not  related  to  the  critical  failure  load.  As  a  result,  large  discrepancies 
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between  the  theory  and  experiment  can  result. 

A  study  on  the  selection  of  appropriate  mathematical  models  for  different 
unidirectional  composite  systems  was  made  [2]  in  the  case  of  static  loading* 

Many  of  the  assumptions  in  [2]  will  also  be  used  in  the  dynamic  problem  treated 
here.  One  of  them  is  the  existence  of  inherent  flaws  or  cracks  which  are  the 
sites  of  failure  initiation. 

Analytical  investigation  of  the  fracture  of  fibrous  composite  materials  sub¬ 
jected  to  impact  loading  has  been  meager  because  the  elastodynamic  stress  analy¬ 
sis  involves  numerous  parameters  and  is  enormously  complex.  This  is  necessitated 
by  the  complex  nature  of  the  dynamic  load  transfer  characteristics  in  composites 
containing  initial  imperfections  such  as  flaws  or  cracks.  The  stress  wave  solu¬ 
tion  is  not  only  time-dependent  but  it  interacts  with  the  material  properties  of 
the  constituents  of  the  composite  and  the  various  geometric  parameters.  The  in¬ 
fluence  of  these  parameters  will  be  analyzed  in  this  impact  study  with  particular 
emphasis  placed  on  determining  the  dynamic  stress  intensity  factors  k^  and 
arising  from  normal  and  shear  loading.  Their  combination  (off-axis  or  angle 
loading)  determines  the  response  to  loading  of  a  more  general  nature  and  reflects 
the  energy  absorption  property  of  the  composite.  Several  examples  of  how  k-j  and 
k2  can  be  combined  to  predict  crack  behavior  in  dynamic  stress  fields  are  found 
in  [3].  The  question  of  whether  there  is  the  need  of  how  to  define  a  dynamic 
fracture  toughness  parameter  differing  from  its  corresponding  static  value  has 
been  the  subject  of  many  past  and  present  debates  within  the  fracture  mechanics 
community.  Thus  far,  no  general  agreement  has  been  achieved. 
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This  report  is  concerned  with  dynamic  fracture  analysis  and,  particularly, 
with  the  development  of  an  analytical  method  for  obtaining  effective  dynamic 
stress  solutions  to  unidirectional  composites  with  cracks  embedded  in  the  ma¬ 
trix.  Other  possible  failure  modes  will  be  dealt  with  in  future  reports.  Ef¬ 
fective  stress  solutions  for  and  k2  are  essential  as  they  are  the  prerequi¬ 
sites  for  formulating  failure  criteria  and  guiding  the  development  of  numerical 
procedures. 


ANGLE  CRACK  UNDER  IMPACT 

Figure  1(a)  considers  a  crack  in  a  layer  of  matrix  material  of  thickness 
2h.  The  composite  is  reinforced  by  unidirectional  fibers  that  are  aligned  par¬ 
allel  with  one  another  and  make  an  angle  with  the  time-dependent  applied  stress 
cT(t).  Without  serious  loss  in  generality,  the  composite  is  assumed  to  be  modeled 
by  a  layer  of  cracked  material  with  elastic  properties  ,  v-j  and  p-|  sandwiched 
in  between  two  dissimilar  media  with  properties  y2,  '^2  ^2’  Ub).  The 

number  of  layers  surrounding  the  cracked  layer  is  reasonably  large  so  that  the 
average  shear  modulus  y2,  Poisson's  ratio  V2  and  mass  density  P2  can  be  used. 

The  basic  two-dimensional  elastodynamic  equations  in  the  theory  of  elastic¬ 
ity  can  be  expressed  in  terms  of  two  scalar  potentials  4)-(x,y,t)  and  Tiij(x,y,t) 

J  J 

where  i,j  =  1,2  with  1  and  2  referring  to  the  cracked  layer  and  the  surrounding 
material,  respectively.  In  terms  of  the  Lamd  coefficients  X.  and  y.,  the  dy- 

J  v) 

namic  stress  components  are 


-4- 


92(j)  32^ 

(Ox)  =  +  2«J  (33^  +  Jjji) 

J 


'“y> .  ' 

-J 


X.  X.+2u. 

(^.)  .  =  /  (I^) 


32^. 

(Txy) ,  =  Uj  (2  ^  ^  -  33P^) 

J 


(1) 


where  V^  =  3^/3x^  +  3^/3y^  and  the  thickness  shear  stresses  are  assumed  to  van¬ 
ish.  The  corresponding  in-plane  displacements  are  given  by 


(2) 


Under  plane  strain,  the  material  elements  are  constrained  in  the  z-direction. 
The  governing  differential  equations  can  then  be  obtained  from  the  equations  of 
moti on : 


1 


3F" 
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The  problem  involves  the  determination  of  the  potentials  (j).(x,y,t)  and  4'-i(x>y»t) 

J  J 

in  equations  (3)  from  the  transient  boundary  conditions  of  the  crack  problem. 

The  analysis  may  be  simplified  considerably  if  the  problem  is  separated  in¬ 
to  two  parts.  The  first  concerns  with  normal  stresses  applied  to  the  crack  such 
that  symmetry  prevails  about  the  x-axis  in  Figure  1(b)  while  the  second  deals 
with  shear  surface  tractions  so  that  the  problem  is  skew-symmetric  with  refer¬ 
ence  to  the  x-axis.  Both  of  these  problems  will  be  presented  separately. 

NORMAL  IMPACT 

Let  the  composite  body  be  initially  at  rest  such  that  the  stresses  are  zero 
everywhere.  Suddenly,  at  t=0,  a  normal  stress  of  magnitude  -cr^  is  applied  to 
the  top  and  bottom  crack  surfaces  in  Figure  1(b)  and  kept  on  the  crack  of  length 
2a  thereafter.  Referring  to  the  set  of  axes  x  and  y  that  are  placed  parallel 
and  normal  to  the  line  crack,  the  following  conditions  are  prescribed: 


(aj  (x,o,t)  =  -a^H(t);  (t  ,,)  (x,o,t)  =  0,  0<x<a;  t>0 
y  1  u  ^y  1 


where  H(t)  is  the  Heaviside  unit  step  function.  The  symmetry  conditions  about 
the  axis  y=0  are  enforced  by  noting 


(u  )  (x,o,t)  =  0;  (t  )  (x,o,t)  =  0,  x>a;  t>0  (6) 

y  1  xy  ^ 
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Perfect  bonding  will  be  assumed  along  the  interfaces  between  material  1  and  ma¬ 
terial  2.  This  requires  the  stresses  and  displacements  to  be  continuous  across 
y  =  ±h.  On  account  of  symmetry,  only  the  upper  half  plane  y>0  need  to  be  con¬ 
sidered,  i.e.. 


(o-,,)  (x,h,t)  =  (a  )  (x,h,t) 

y  1  y  2 

(7) 

(^xy)^  (X,h,t)  =  C'r^y)^C5^>h»t) 


and  for  the  stresses  and 


(u^)  (x,h,t)  =  (u  )  (x,h,t) 

X  1  X  2 

(8) 

(u  )  (x,h,t)  =  (u  )  (x,h,t) 

y  1  y  2 

for  the  displacements. 

VuuoUi  lY\;t2.QnjaJL  zquatloyL6.  It  is  convenient  at  this  point  to  apply  the 
Laplace  transform  to  the  time  variable  t  which  corresponds  to  p  in  the  trans¬ 
formed  plane.  Consider  the  standard  Laplace  transform  on  f(t): 

f*(p)  =  /  i'(T)  e"*^^dt  (9) 

0 

whose  inversion  is 


f(t)  -  i  /  f*(p)  e^^dp 
Br 


(10) 
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in  which  Br  stands  for  the  Bromwich  path  of  integration.  The  application  of 
equation  (9)  to  equations  (3)  yields 


J  j 


J  C2j  j 


Again,  the  condition  of  symmetry  requires  only  the  consideration  of  x  and  y  in 
the  first  quadrant.  The  Fourier  cosine  and  sine  transforms  defined  by 


f^(s)  =  /  f{x)  cos(sx)dx 


=  f  /  f^(s)  cos(sx)ds 
0 


f^(s)  =  /  fCx)  sin(sx)dx 


=  f  /  f^(s)  sin(sx)ds 
0 


are  now  applied  to  the  space  variable  x.  This  simplifies  equations  (3)  to  a  set 
of  ordinary  differential  equations  which  can  be  solved  giving 


<|)i(x,y,p)  =  f  7  [A^^^s,p)e  +  A^^^s,p)e’^^^'^]  cos(sx)ds 
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(14) 


Pf(x,y,p)  =  f  /  [B^^^(s,p)e  (s,p)e^^^^]  sin(sx)d; 


for  the  cracked  matrix  and 


(i)5Cx,y,p)  =  I  /  C^^^(s,p)e  cos(sx)ds 


4^|(s,y,p)  =  I  /  C^^^(s,p)e  sin(sx)ds 


for  the  averaged  fiber-matrix  material.  In  equations  (14)  and  (15),  the  quanti¬ 
ties  Yij  and  Y2j  ^re  given  by 


n2  1/2  n2  ■'/2 

'  1='  %>  ’  ^2j  =  %’ 


The  functions  in  equations  (14)  and  (15)  are  deter¬ 
mined  from  the  transient  boundary  conditions.  To  this  end,  equations  (5)  and 
(6)  will  be  written  in  the  Laplace  transform  plane: 


(a*)  (x,o,p)  =  -  Cx*  )  (x,o,p)  =  0,  0<x<a 
y  1  p  AY  1 


(u*)  (x,o,p)  =  0;  (t*  )  (x,o,p)  =  0,  x^a 
y  1  xy  ^ 


In  the  same  way,  equations  (7)  become 


-9- 


(19) 


(a*)  Cx,h,p)  =  Cct*)  (x,h,p) 

y  1  y  2 


(t*  )  Cx,h,p)  =  (t*  )  (x,h,p) 
xy  ^  xy  2 


and  equations  (8)  take  the  forms 


(u*)  (x,h,p)  =  (u*)  (x,h,p) 

X  1  X  2 

(20) 


(u*)  (x,h,p)  =  (u*)  (x,h,p) 

y  1  y  2 


The  stresses  and  displacements  in  equations  (1)  and  (2)  may  also  be  trans¬ 
formed  into  the  Laplace  transform  plane.  Without  going  into  details,  the  appro¬ 
priate  Laplace  transform  of  the  stress  and  displacement  expressions  in  equations 
(17)  to  (20)  may  be  used  to  satisfy  all  of  the  necessary  boundary,  symmetry  and 
continuity  conditions.  This  leads  to  the  following  set  of  dual  integral  equa¬ 


tions; 


00 

/  A(s,p)  cos(sx)ds  =  0,  x^a 
0 

(21) 

Tra 

/  sFj(s,p)  A(s,p)  cos(sx)ds  =  -  4p^p(Y,Kp'» 


in  which  Fj(s,p)  stands  for  the  known  function 
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'  s(1-Kf)A^  -  6<^>e  ] 

^  [1  (s^n|,)“  ^  (22) 

while  the  quantities  k-j  and  Aj  are  defined  as 

1  /2 

K]  =  (c2i/Ci 1 ) 

(23) 

+  B<') 

such  that  6*^'.  6*^^ — .  6  ^  are  given  by 

gd)  ,  („(3),(6)  .  „C2)„(7))/a^;  6(2)  =  (at'Hc.'®)  - 

(24) 

6(3)  .  (a(')c.(3)  .  c(3)a(5))M„;  B<‘'>  =  (o.'^>a(8)  -  .('‘)c.(5))M^ 


where  A^  is 

A  =  (25) 

0 

The  quantities  in  equations  (25)  are  complicated  functions 

of  s,  p  and  the  material  constants.  They  are  given  by  equations  (I.l)  in  Appendix 


The  problem  is  now  reduced  to  finding  the  single  unknown  A(s,p)  governed 
by  equations  (21).  Once  A(s,p)  is  known,  the  functions  A^^^,  A^^^, — , 
that  are  required  in  equations  (14)  and  (15)  for  the  Laplace  transform  of  the 
potentials  (})T(x,y,p)  and  4't(x,y,p)  can  be  obtained  from  equations  (1.2)  outlined 

J  J 

in  Appendix  I.  What  remains  is  the  determination  of  a  solution  for  the  dual  in¬ 
tegral  equations  (21).  This  will  be  accomplished  with  the  help  of  a  method  by 
Copson  [5]  which  has  been  used  by  Chen  and  Sih  [6]  for  solving  dynamic  crack 
problems  involving  single-phase  homogeneous  materials.  Following  the  details 
in  [5,6],  it  can  be  shown  that 

TTo  a^  1  _ 

=  =  -  4y,p(1-.p  I  ''5  ♦!(«•'>>  (26) 

is  a  solution  of  equations  (21)  with  being  the  zero  order  Bessel  function  of 
the  first  kind.  The  function  $|(C,p)  is  calculated  numerically  from  a  Fredholm 
integral  equation  of  the  second  kind: 

1 

KT(?,n,p)dn  =  vf  ■  (27) 

1  Q  1  1 

whose  kernel 


Kl(5.n.p)  =  /  s[Fj(-|,  p)  -  1]  Jq(sC)  JQ(sn)ds  (28) 

0 

is  symmetric  in  ^  and  n- 

Mode  I  dynamic,  ipiui  Inttn^iity  ^acto^.  The  transmission  of  the  time-de¬ 
pendent  load  to  the  vicinity  of  the  crack  tip  can  be  best  described  by  the  in¬ 
tensification  of  the  local  stresses.  A  quantity  that  has  been  used  widely  in 
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the  static  theory  of  fracture  mechanics  is  the  "stress  intensity  factor"  which 
can  be  extracted  from  the  asymptotic  expansions  of  the  stresses  near  the  crack 
tip.  Referring  to  Figure  2,  let  r-j  and  6i  be  a  set  of  local  polar  coordinates 
measured  from  the  right  hand  crack  tip  located  at  x=a  and  y=0  in  the  matrix 
material,  the  singular  character  of  the  dynamic  stresses  is  described  only  by 
the  space  variables  and  hence  can  be  more  easily  determined  in  the  Laplace 
transform  domain.  This  observation  was  first  made  by  Sih,  Ravera  and  Embley 
[7].  Following  their  procedure,  the  local  stresses  in  terms  of  r^  and  are 
found: 


k*(p)  e,  01  30, 

(ot^)  (r,,0,,p)  =  cos  ^  (1  -  sin  j-  sin  -y)  +  0(r|) 

X  1  1  I  /2r^ 


k?(p)  0,  01  30, 

(a*)  (r,,0,,p)  =  -j - cos  y  (1  +  sin  ^  sin  — )  +  0(rp 

y  1  1  1  ^  ^ 

(a*)  (r,,0,,p)  =  4iir2v,  cos  y  +  0(rp 
z  1  I  I 

k?(p)  0,  01  30, 

Ct*  )  (rp0^,p)  =  sin  y  cos  y  cos  -y  +  0(rp 


(29) 


Only  the  dynamic  stress  intensity  factor,  k|(p),  in  equations  (29)  need  to  be 
inverted  to  real  time  t: 


kf(p) 


$|(l,p) 

P 


(30) 


where  the  function  $|(l,p)  is  found  from  $|(5,p)  by  letting  the  nondimensional 
parameter  C=1  representing  the  crack  tip  location.  The  functional  dependence 
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of  the  stresses  in  r-j  and  S-j  as  shown  by  equations  (29)  reveals  that  the  dynamic 
stresses  also  possess  the  inverse  square  root  singularity  in  terms  of  r-j  and 
that  the  angular  distribution  in  9,  is  the  same  as  the  case  for  static  loading. 

Applying  the  Laplace  inversion  formula  in  equation  (10)  to  (30)  renders 
the  factor  k^(t)  as  a  function  of  time,  i.e.. 


k^(t)  = 


a^/a  _  $|(l,p) 


2iTi 


/ 

Br 


e^^dp 


(31) 


It  is  apparent  that  $*(l,p)  must  be  first  known  before  the  integration  of  equa¬ 
tion  (31)  can  be  performed.  Refer  to  Appendix  II  for  a  detailed  account  of  the 
procedure  used  for  evaluating  equation  (31).  Three  different  sets  of  <I>|(l,p) 
values  are  plotted  against  the  dimensionless  Laplace  transform  wave  number  C2i/pa. 
They  are  given  in  Figures  3  to  5  for  Pi  =  P2  and  v-j  =  ^2  “  while  the  ratios 
a/h  and  P2/fi-|  varied.  In  general,  all  the  curves  tend  to  rise  quickly  and 
then  flatten  out.  It  would  be  more  meaningful  to  discuss  the  influence  of  a/h 
and  v-2fv^-\  the  stress  intensity  factor  ki(t). 

Figures  6  to  8  display  the  normalized  stress  intensity  factor  k^(t)/a^/a 
as  a  function  of  C2it/a.  In  Figure  6,  the  crack  length  to  layer  thickness  ratio, 
a/h,  is  fixed  at  unity  while  the  shear  moduli  ratio,  V2^V‘'\  increased  from  0.1 
to  10.0  as  indicated.  The  k-j  (t)  factor  is  oscillatory  in  nature  reaching  a  peak 
and  then  decreases  in  magnitude  as  time  increases.  The  oscillation  is  more  pro¬ 
nounced  wheii  the  shear  modulus  of  the  cracked  material  is  greater  than  that  of 
the  surrounding  material,  i.e.,  y-j  >  y2.  The  values  of  k-j  (t)  decrease  below 
those  of  the  corresponding  homogeneous  case,  y^  =  y2,  solved  previously  by  Chen 
and  Sih  [6]  when  y-j  <  y2.  The  influence  of  a/h  on  ki(t)  is  exhibited  in  Figures 
7  and  8  for  the  two  cases  of  y2/yi  =0.1  and  10.0,  respectively.  For  y2/yi  =  0.1 
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in  Figure  7,  a  decrease  in  a/h  tends  to  lower  the  stress  intensity  factor.  Ob¬ 
served  also  is  a  small  step  in  the  curve  for  a/h  =  2.0  and  small  time  t.  This 
corresponds  to  the  reflection  of  elastic  waves  from  the  material  interface. 

The  size  and  time  scale  are  such  that  this  effect  showed  up  quantitatively  in 
the  graph  while  the  same  effect  was  not  noticeable  in  the  other  curves.  For 
the  smaller  ratios  of  a/h  such  as  0.5  and  1.0,  the  crack  tips  are  further  away 
from  the  interface  and  the  influence  of  the  reflected  waves  are  not  as  pro¬ 
nounced.  The  opposite  trend  is  observed  in  Figure  8  for  V2/P1  =  10.0.  When  the 
outer  material  is  more  rigid  than  that  of  the  center  layer,  k.|(t)  tends  to  in¬ 
crease  in  magnitude  as  a/h  is  decreased.  Again,  a  distinct  step  in  the  curve 
for  a/h  =  2.0  is  seen  for  small  time  t.  As  time  increases,  all  of  the  results 
here  reduce  to  the  corresponding  static  solutions  of  Hilton  and  Sih  [8]. 

Gm2Aal  loading.  If  the  normal  stress  applied  to  the  crack  surface  is  not 
constant  in  magnitude  but  may  vary  as  a  function  of  x,  then  the  dynamic  stress 
intensity  factor  can  be  obtained  by  adding  a  sequence  of  solutions  corresponding 
to  step  loadings  with  different  stress  levels  Oq,  etc.  In  other  words,  the 
general  loading  a(t)  may  be  considered  as  the  sum: 

(j(t)  -  OgH(tg)  +  (JiH(ti)  +  +  ...  (32) 

This  is  illustrated  graphically  in  Figure  9.  From  equations  (31)  and  (32),  the 
factor  k.|(t)  that  corresponds  to  a{t)  may  be  written  down  inmediately  as  follows 

,  H(1>P)  nt 

*^1^^^  ^  Tin  '*■  ^l^^^l^  p  ® 

Equation  (33)  may  be  used  to  derive  k.|(t)  for  any  time-dependent  normal  surface 

tractions  which  in  turn  can  also  simulate  any  loadings  that  are  applied  at  dis- 
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tances  away  from  the  crack  by  means  of  the  principle  of  superposition. 

SHEAR  IMPACT 

Suppose  that  the  crack  in  Figure  1(b)  is  now  sheared  suddenly  by  a  pair  of 
shear  stresses  of  magnitude  such  that  the  upper  and  lower  crack  surfaces 
move  in  the  opposite  direction.  This' creates  a  deformation  field  that  is  skew- 
symmetric  with  respect  to  the  y=0  plane.  Following  the  footstep  laid  out  in  the 
previous  example  on  normal  impact,  the  Laplace  transform  of  the  transient  bound¬ 
ary  conditions  on  the  x-axis  inside  the  crack  are 


(t*  )  (x,o,p) 
xy  ^ 


(x,o,p)  =  0,  0<x<a 


and  the  skew-symmetric  conditions  outside  the  crack  are  given  by 


(34) 


(ut)  (x,0,p)  =  0;  (ct*J  (x,o,p)  =  0,  x>a 


(35) 


1 


1 


Continuity  of  the  stresses  across  y=h  is  expressed  by 


(a*)  (x,h,p)  =  (a*)  (x,h,p) 

y  1  y  2 


(36) 


(t*  )  (x,h,p)  =  (t*  )  (x,h,p) 


xy 


1 


xy'2' 


while  the  displacements  are  also  required  to  be  continuous,  i.e.. 


(u*)  (x,h,p)  =  (u*)  (x,h,p) 

X  1  X  2 


(u*)  (x,h,p)  =  (u*)  (x,h,p) 


(37) 


1 


y' 


-16- 


Tyvtzgfutt  H.tpfiz^^nt(xtlQYi&.  Under  the  above  considerations,  the  following 
wave  potentials  <j)^(x,y,p)  and  M(x,y,p)  are  selected: 

J  J 


4>f(x,y,p)  =  f  /  [A^^^(s,p)e  sin(sx)ds 

0 

(38) 

iiJf(x,y,p)  =  I  /  [B^^^(s,p)e  +  B^^^(s,p)e^^^^]  cos(sx)ds 


for  the  cracked  layer  and 


o  “  / 1  \  “  Y  i  py 

<j)f(x,y,p)  =  f  /  C^‘Hs,p)e  sin(sx)ds 
0 

(39) 

4'2(x>y>P)  =  I  /  C^^^(s,p)e  cos(sx)ds 


for  the  outside  material. 

Equations  (38)  and  (39)  may  now  be  substituted  into  the  Laplace  transform 
of  the  stresses  and  displacements  in  equations  (1)  and  (2).  Making  use  of  the 
conditions  in  equations  (34)  to  (37),  the  solution  can  be  expressed  in  terms  of 
the  functions  A^^^,  which  are  related  to  a  single  unknown  B(s,p) 

as  shown  by  equations  (III.l)  in  Appendix  III.  The  function  B(s,p)  is  governed 
by  the  system  of  dual  integral  equations 


oo 

/  B(s,p)  cos(sx)ds  =  0,  x^a 
0 


OO  TTT 

/  sFjj(s,p)  B(s,p)  cos(sx)ds  =  x<a 


(40) 
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The  function  FjjCs,p)  is  related  to  Fj(s,p)  in  equation  (22)  as 


Fll(s,p)  =  Fj(s,p) 


where 


^II  2^^  '21' 


The  other  parameters  such  as  Ky  Aj,  3^^^,  etc.,  are  the  same  as  those  de¬ 

fined  earlier  for  the  case  of  normal  impact. 

A  solution  to  equations  (40)  is  again  found  by  application  of  the  Copson's 
method  [5]: 


TTT  a^  1  _ 

■  4p,p(1-<f)  I  ''5  *!l(5'P) 


provided  that  $|j(?,p)  satisfies  a  Fredholm  integral  equation  of  the  second  kii 


'*’*l(?.p)  +  /  '^’*i(n,p)  Kjj(?,n,p)dri  =  i/f 


whose  kernel  Kjj(^,n,p)  takes  the  form 


Kjj(?,ri.p)  =  !  s[Fjj(|-,  p)  -  1]  Jg(s?)  jQ(sn)ds 

0 


Made  II  dynamic.  Intzn^lty  ^aato^.  As  in  the  case  of  Mode  I,  the 

asymptotic  expressions  of  the  dynamic  stresses  in  the  Laplace  transform  plane 
are  first  obtained  in  terms  of  r-j  and  O-j  defined  in  Figure  2.  The  results  .are 
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ko(p)  01  01  36, 

(a*)  (r,  ,9,  ,p)  =  -  — - sin  ^  (2  +  cos  ^  cos  -^)  +  0(rp 

X  ^  I  I  ^  ^ 

l^2^P)  01  ®i  36, 

(a*)  (r,  ,9,  ,p)  =  - - sin  ^  cos  ^  cos  -s-  +  0(r, ) 

y  1  I  I  ^  6  c  c  \ 

(46) 

^2^^^  0] 

(ct*)  (r,,9,.p)  =  -  -:^2v,  sin  0(rp 
2  1  1  I  ^  I  ^  I 

koCp)  01  ®i  39, 

(t*  )  (r,  ,6,  ,p)  =  — - sin  ^  cos  ^  cos  -y"  +  0(r, ) 

1  '  '  /2r^  2  1 

with  k5(p)  being  the  only  quantity  that  depends  on  time  through  the  parameter  p: 


^frd.p) 

k|(p)  =-^^p - T 


Equation  (10)  is  then  applied  to  invert  the  Laplace  transform  of  the  stress  in¬ 
tensity  factor  in  equation  (47).  This  gives 


k2(t)  = 


in  which  $|j(1jP)  is  computed  numerically  from  equation  (44). 

Figures  10  to  12  display  the  values  of  $^j(l,p)  as  a  function  of  the  normal¬ 
ized  quantity  C2i/pa  for  various  values  of  a/h  and  ^2^^!  while  v.|  =  V2  =  0.29  and 
p.|  =  p2  are  used  for  all  cases.  With  a  knowledge  of  ^>|j(l.p)j  the  integral  in 
equation  (48)  may  be  evaluated  by  a  procedure  outlined  in  Appendix  II.  In  gen¬ 
eral,  k2(t)  increases  with  time  reaching  a  maximum  and  then  decreases  to  the 
static  value  for  sufficiently  large  time.  The  trend  is  very  similar  to  k-jCt) 
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for  the  case  of  normal  impact  in  that  a  higher  value  of  k2(.t)  is  obtained  when 
the  modulus  of  the  surrounding  material  is  less  than  that  of  the  cracked  layer, 
i.e.,  y2/^l  Comparing  the  results  in  Figures  6  and  12,  it  is  seen  that  for 

^2^^1  impact  yields  a  higher  crack  tip  stress  intensity  factor  than 

shear  impact,  i.e.,  k.|(t)  >  k2(t).  The  opposite  is  observed  when  v-2/v-^  ^  i.e., 

k2(t)  >  k-|(t).  The  curves  in  Figures  14  and  15  for  k2(t)  also  show  the  absence 
of  a  small  fluctuation  for  small  time  which  was  present  in  Figures  7  and  8  for 
k.|(t).  This  is  because  the  influence  of  the  reflected  incident  shear  wave  from 
the  interface  is  considerably  weaker  even  for  the  ratio  of  a/h  =  2.0. 

CONCLUSION 

As  composite  materials  are  currently  being  applied  to  major  primary  structure 
designs,  it  is  necessary  to  have  an  in-depth  understanding  of  the  mechanical  be¬ 
havior  of  these  materials,  particularly  with  reference  to  the  allowable  applied 
load  both  statically  and  dynamically.  This  investigation  is  concerned  with  the 
dynamic  stress  distribution  around  a  crack  embedded  in  the  matrix  of  a  unidirec¬ 
tional  composite.  The  time-dependent  loading  can  be  of  a  general  nature  applied 
in  an  arbitrarily  direction  with  reference  to  the  crack  plane.  For  those  com¬ 
posites  which  fail  predominantly  by  matrix  cracking  under  impact,  the  present  re¬ 
sults  can  be  used  effectively  for  determining  the  ability  of  the  composite  to 
absorb  energy  and  to  withstand  load  prior  to  total  destruction. 

The  other  modes  of  failure  such  as  fiber  breaking  and/or  debonding  of  fibers 
from  matrix  are  not  treated  but  may  be  significant  in  other  composite  systems. 

The  redistribution  of  dynamic  stresses  in  these  cases  may  also  be  analyzed  such 
that  their  individual  contribution  can  be  assessed  quantitatively.  These  cases 
will  be  left  for  future  investigations. 
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APPENDIX  I:  EXPRESSIONS  FOR  AND  A^^^(s,p),  — ,  C^^^(s,p) 

IN  NORMAL  LOADING 

This  section  gives  the  expressions  for  in  equations 

(25)  in  terms  of  the  variables  s,  p  and  the  material  constants 


(11  ^2  /  p2  w  ^2r^22  N-, 

=  -s[{l  -  — )t9i  -  —  V2^)(s^-v 


^1  ^^22  ^  "^12^22^ 


(2)  .  .rn 

a  -  s[(l  y^)Y21  ^ 


J3)  =  1  (sa^2  )  (^7-^)] 

^  2  ^2r  y-j  2c^2 


12'22 


,(4)  .  1  .^yu-<22 

a  -  9  ^2r  y-  U^-v. -v-- 


2c|2  s‘-y-|2T22 


)] 


cC5)  ,  .  1  (,2+^2  w  rs2  +  ^  ("'^1,2121 

a  2  ^^21  ^  ^  y^  2^  4"-Yi2Y22 


)] 


(I.l) 


-ns,  1.5 1...  *,S5gi,i 


22  ^  ■^12‘^22 


•  «■  ■  -  ?  '4’'i5Sfe" 


(8)  _  .rn  ^2> 

'  --'r"  -i:7”ll  -^7'2i|2''i^T,2Y22 


'2  ,  p2  ■^ll'^^12  ^^ 
a-  '  =  -s[(l  -  — )Yn  -  —  (9^)(<.2_^:_---)] 


in  which  y,-,-  is  given  by  equations  (16). 

•  J 
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The  functions  A  ^  are  related  to  the  single  function  A(s,p) 


as  follows: 


A<'>(s.p)  = 
A<^'(s,p)  = 


e''"2'')-sy 


At  ‘■2  ^  ^  ^  ~ 


rsY  e 


■(Yll+Y2,)h  ,  1 


-2y 


n 


-(Y-„tY2,)h 
e  J 

'  (b'" 


+  e 


(3) 


h 

)] 


B<l)(s.p)  =  S<'>  a(1)(s.p)  t  ^(2)  A<2)(s.p) 

b(2)(s,p)  =  e(3)  a(i)(3.p)  ,  ^(4)  J^ir^si)'’ 

“  F-7i1y32  [<''-flll'22)  ®  a('>(s.p) 

+  (4^+YiiY22)  e  "  A^2)(5_p)  +  sCyji-Yzj)  ®  B^''(s,p) 

-  s(Y21^22>  b'^'Cs.p)] 

-  4(y,i+Y22)  A*^'(s,p)  +  (4“-Ti2Y2i'  ®  B^'^s.p) 

t  (s"ni2Y2l)  b'^>(s.p)] 


(1.2) 
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APPENDIX  II:  METHOD  FOR  EVALUATING  THE  DYNAMIC  STRESS 
INTENSITY  FACTOR  EQUATION  (31) 

The  integral  in  equation  (31)  is  basically  of  the  form 


g(t)  =  I 


Br 


dp 


(II. 1) 


The  Bromwich  path,  Br,  consists  of  an  infinite  line  parallel  to  and  to  the  right 
of  the  imaginary  axis  in  the  complex  p-plane.  The  function  f*(l,p)  is  considered 
to  be  known  for  discrete  values  of  p.  There  are  a  number  of  available  methods 
for  finding  g(t)  as  a  process  in  the  Laplace  inverse  transform.  The  method 
adopted  here  can  be  found  in  [9,10j. 

The  integral  f*(l,p)/p  in  equation  (II. 1)  is  first  evaluated  at  the  points 


p  =  (l+n)(S,  n  =  0,1 ,2, — 


(II. 2) 


in  which  5  is  a  real  and  positive  number.  According  to  equations  (9)  and  (10), 
f*(l»p)/p  may  be  written  as' 

111^=  /  g(t)e'P^dt  (II. 3) 

P  0 

The  above  infinite  integral  is  now  transformed  to  a  finite  integral  on  the  in¬ 
terval  [-1,1]  by  making  the  substitutions 


*f*(l,p)  stands  for  $|(l,p)  in  normal  impact  and  3>|j(l,p)  in  shear  impact  and 

they  are  calculated  from  the  Fredholm  integral  equations  of  the  second  kind, 
namely  equations  (27)  and  (44). 


-23- 


1 


(II. 4) 


X  =  2e" 


S(x)  =  g[t(x)]  =  g[-  1  log(^)] 


:ii.5) 


Therefore,  equation  (I I. 3)  becomes 


=  -irr  I  O+x)'^  G(x)dx 


2n+1 


(II. 6) 


in  which  G(x)  can  be  expanded  in  series  form  consisting  of  Legendre  polynomials 
P^(x)  which  are  orthogonal  on  the  interval  [-1,1],  i.e.. 


G(x)  =  I  C.  P.(x) 
i=0  ^  ^ 


(II. 7) 


Similarly,  the  function  (l+x)*^  in  equation  (II. 6)  may  also  be  expanded  in*  the 


(l+x)"  =  I  D.  P.(x) 
i=0  ^  ^ 


(II. 8) 


such  that 


n  -  on/o,.  n(n-l )— [n-(i-l )! 
°i  "  2  (n+l)(n+2)---  n+iil 


(II. 9) 


Putting  equations  (II. 7)  and  (II. 8)  into  (II. 6)  and  applying  the  orthogonality 
conditions  for  the  Legendre  polynomials,  the  following  sum  is  established: 
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(II. 10) 


•Thus  the  coefficients  C^.  may  be  found  with  given  by 


Cq  =  f*Cl,6) 


(II. 11) 


For  a  finite  number  of  N  coefficients,  a  partial  sum  for  G(x)  in  (II. 7)  is  ob¬ 
tained  and  an  approximate  evaluation  of  g(t)  can  be  made  since  from  equation 
(II. 5) 


g(t)  =  I  C.P.  LZe'^^  -  1] 
i=0  ^  ^ 


(11.12) 


The  parameter  5  is  chosen  such  that  g(t)  is  best  described  for  the  range  of  t 
considered. 
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APPENDIX  III:  EXPRESSIONS  FOR  A^^^s,p),  — ,  C^^^(s,p) 
IN  SHEAR  LOADING 


In  the  skew- symmetric  problem,  the  unknown  functions  in  equations  (38)  and 
(39)  can  also  be  expressed  in  terms  of  a  single  unknown  B(s,p)  in  accordance  with 
the  following  relationships: 


A^^^s.p) 


=  .  B(s,p) 
^II 


[SY21 


(2)  _ 


)  +  \  (s^V 


21 


)e 


-(Yii+T2i)h 


] 


A<2)(s,p)=^[sY2,e''""'  (6'’>  -  6<3V‘’21") 


-2Yoih. 


loo" 

+  \  (s^+Y2-i  )e  ] 


B^^^(s,p)  = 


B(')e 


(Yii-Y2i)h 


21 


)h 


.(2) 


(s,p) 


B<2)(s,p)  =  -  a<1)(s.p)  -  a'^'Is.p) 

/■|\  o^^2*^  “Yiih  /tn 

^  "  ^''"Xi2^22  ^^^^"^11^22)®  ^ 


+  (s^+YiiY22)e  A^^^(s,p)  -  s(y2-]-Y22)®  B^^^(s,p) 

+  s(Y2T+Y22)e^^^  B^^^(s,p)] 
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(III.l) 


+  s(Y^2'^n)®^^^  A^^^(s,p)  +  (s^-Y2iY-i2^®  B^^^(s,p) 

+  {s^+Y2iY]2)e^^^  B^^^(s,p)] 

where  Ajj  is  given  by  equation  (42). 
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Figure  2.  Stress  element  near  crack  in  matrix  of  fiber-reinforced 
composite 
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k|(t)/cr, 


2.4 


a/h  =  1.0 
V\  =  1/2  =  0.29 

P\  =Pz 


Figure  6.  Dynamic  stress  intensity  factor  k, (t)  versus  time  for 
a/h  =  1.0 
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a/h  =  2.0 


M2//^i  =  0.I 
i/|  =  z/2  -  0.29 


C  2  I  t/Q 


Figure  7.  Dynamic  stress  intensity  factor  ki(t)  versus  time  for 


Figure  8.  Dynamic  stress  intensity  factor  ki(t)  versus  time  for 
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C(i.p) 


cr 


Time  t 


Figure  9.  Applied  stress  as  a  general  function  of  time 


Figure  10.  Variations  of  $jj(l,p)  with  C2i/pa 
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a/h-  2.0 


4^ 

Figure  11.  Variations  of  «jj(l»p)  with  C2-j/pa 


a/h  =  0.5 


*  10.0 
i/|  s  ^2  "  0.29 


P\  ^Pz 

I  I _ I - 1 - L. 

1.0  2.0  3.0  4.0  5.C 

Cj/pa 

Figure  12.  Variations  of  $jj(l,p)  with  C2-|/pa 
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k2(t)/T 


Figure  14.  Dynamic  stress  intensity  factor  k2(t)  versus  time  for 
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COMPUTER  PROGRAMS 


Nonmal  ■mpaat. 


3 

3 

I 

3 

3 

'  I 

.  5 

20 


20 

22 

23 

24 

31 

33 

34 

35 
37 
41 

46 

47 
54 

56 

57 
61 
64 
64 
72 
72 

74 

75 

104 

106 

110 

111 

114 

123 

123 

134 

141 

144 

147 

151 

160 

160 

163 

165 

167 

171 

205 

205 

207 

212 


program  B|TA<roPUT,OUTPUT,PUNCH,PLOI,TAPE 
REAL  N0N(4)  ,F(4*4»1)  ,6(4»4J  ,uiH#  tri  at 

iJt  e^tj5?.'STA<50, 

equivalence’ (NON, B) 

CCMMON  t<l»  K2fK3jK4  ^ 

COMMON/ AUX/HtP»P'<l»P*^2»  0MU»  X,t 

LPtDfO.O 

--  REA0^2» K1»K2,K3»K4 

♦  -°oroer^cf  system  of  equations 

:  ;  NO?  OF  glillNCT  KERNELS 

:  ;  iJS:  OF  DATA  sItS  TO  BE  EVALUATED 

♦  SET  UP  DATA  POINTS 

AK=K3 

00  5  N=1»K3 
AN=N 

•  i£T'’up''iNTEGRATION  MATRIX 

M-K3-2 

N=K3-1 

A=1^/(3.’A) 

DO  10  K=2*M»2 
10  0(K)=2.*A 

□0  15  K=1,N,2 
15  DCK)=4.*A 

*  CALCULAli  NONHOMOGENEOUS  TERMS 

RHS=1.0  ^ 

00  22  I=1»K2 
PRINT  9 
''  FORMAT (IHl) 

READ  61xBMU 
FORMATCFIO. 5) 

00  999  n  =  l»K4 

o=-  non  (N)  =  Rh|*SQRT  (PT(N)  ) 

'='“-grL£^So!SI?aV 

00  20  N=1»K3 
00  20  M=1,K3 
IF(M-N)  25»3a*30 
FCMtN  tl) =F (N  tM» I) 

30  F(hTn,I>=FU(1,PT(MI,PT(N)) 

continue 

liXl  £?Sp^r55?0?’“K3. 

PRINT°6*iFT|L)  ,NCN|L) 

FORMAT  C5X»F8*4»F15.6 1 

fg^Iiflf  =  N0N(K3. 

OTA  (II+1>=P 

PUNCH^6li  (DTAdX)  tLPCIX)»IX=l»19J 

g2TiS«i55-.iiA.LP. 

22  CONTINUE 
END 


99= 


61 


35 


25 


20 


40 


999 

66 


6 

6 

10 

12 

13 

14 
14 
26 
35 


£»Sul5K"^?p5l?PKE.eMU,x.Y 

si'i5?L?ro:EiHo 

45  SIMP=0.0 
RETURN  ^ 

50  CONTINUE 

SA=Z(I»  A) (I»B) 

SB=Z(I, $+2.*0cL) 

SC  =  Z(I»A  +  DiiL)+Z( 


ItA+3.*CEL) 
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PLOT) 


( 

53 

Sl=(0£L/3. <SA+2.*SB+4.’SC) 

61 

IFCSl.EC.O. 0)  GO  TO  45 

62 

K=8 

63 

35 

sa=SB+sc 

-  65 

DEL=0.5*0EL 

67 

SC=Z(I, A*D£L) 

- 75 

_  .  - - 

J=K-1 

77 

DO  5  N=3,J,2 

,  100 

AN=N 

101 

5 

SC=SC+Z (I,A+AN»OEL) 

113 

_ _ - 

S2=(0EL/3,>*(SA*2.*SB+4.»SC) 

— 

122 

DIF=A8S  C(S2-S1)  /SI) 

125 

ER=0. 01 

127 

IF(DIF-ER) 30,25,25 

t 

131 

30 

SIMP=S2 

... 

133 

RETURN 

133 

25 

K=2^K 

134 

S1=S2 

136 

IF  <K-2048)  35,35,40 

14  0 

40 

PRINT  42,I,A,E  ^  ^  , 

152 

42 

F0RMAT(5X,*  INT.  DOES  NOT  CONVERGE  »,I3,2Fg.4} 

152 

PRINT  60,X,Y 

162 

60 

FORMAT<2F10.5) 

—  -• 

162 

DO  70  J=l,10 

166 

DIP  =  J 

167 

DIP=DIP/10. 

171 

H=Z(I,DIF) 

175 

PRINT  60, W 

202 

70 

CONTINUE 

206 

CALL  EXIT 

207 

END 

SUBROUTINE  CHANGE C F, G , D , I > 

7 

REAL  F(4,4,l),Gt4,4) ,0(4) 

7 

COMMON  K1,K2,K3,K4 

7 

DO  10  N=1,K3 

10 

DO  10  M=1,K3 

11 

G(M,N)  =F( M,N,I> ’DCN) 

24 

10 

CONTINUE 

30 

DO  20  N=1,K3 

31 

20 

G(N,N)=G(N,N)+1.0 

40 

RETURN 

41 

END 

SUBROUTINE  LINEQ  ( A  ,E  ,  T ,  N) 

REAL  A(N,N) , 6{N) ,TtN) 

7 

7 

DO  5  I=2,N 

10 

5 

A  (I,1)  =  A  CI,1)/A  {1,1) 

17 

DO  10  K=2,N 

20 

M=K-1 

22 

DO  15  1=1, N 

23 

15 

T(I)=A(I,K) 

33 

DO  20  J=1,M 

34 

A{  J,K)  =T(  J) 

41 

J1=J+1 

43 

DO  20  I=J1,N 

T(I)=T{I)-A{I,J)*A(J,K) 

44 

55 

20 

CONTINUE 

61 

A(K,K)  =  T{K) 

65 

IFCK.EQ.N)  GO  TO  10 

66 

M=K  +  1 

70 

DO  25  I=M,N 

71 

105 

25 

10 

A(I,K)  =  T{I)  /A(K,K) 

CONTINUE 

♦  BACK  SUeSTITUTc 

110 

DO  30  1=1, N 

111 

T(I)=8(I) 

114 

M=I  +  1 

116 

IF(M.GT.N)  GO  TO  30 

121 

DO  30  J=M,N 

122 

B( J)=8{ J) -A { J,I ) ’TCI) 

132 

30 

CONTINUE 

136 

DO  35  I=1,N 
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\ 


137 

K=N4-1-I 

x41 

B{KI=T(K)/A  (K,K) 

146 

K1=K-1 

150 

IF(K1.£Q.0)  GO  TO  35 

151 

00  35  J1=1»K1 

152 

J=K-J1 

154 

T{J)=T( JV-A(J,K)*B(K) 

162 

35 

CONTINUE 

167 

RETURN 

167 

.  .  _ 

END 

FUNCTION  BESJO(A) 

3 

IF<A-3.>5,5,10 

5 

5 

asA*A/9. 

7 

H=l.-2.  2A99997*B 

12 

z=B*e 

13 

H=W+1.265620fi*Z 

15 

Z=Z*B 

16 

W=H-.3163866*Z 

20 

Z=Z*B 

22 

W=W+.0444479»Z 

24 

Z=Z*B 

25 

W=W-.0035444*Z 

27 

Z=2*B 

31 

8ESJO=M+.00021*Z 

34 

RETURN 

34 

10 

e=3./A 

36 

W=. 79788456-. 00000077*6 

41 

\/=A-,78539816-.  041 66397*6 

44 

Z=B*B 

45 

W=W-.0055274*Z 

47 

V=V-.000G3954*Z 

52 

Z=Z*B 

53 

WsW-.  000C95l2*Z 

55 

V=V  +  . 0026257  3*Z 

57 

Z=Z*B 

61 

W=W+.00137237*Z 

63 

V=\I-.  00054125*Z 

65 

Z=Z*6 

67 

W=W-.00072305*Z 

71 

V=V-.00029333*Z 

73 

Z=Z*8 

75 

H=W+. 00014476*Z 

77 

V=V+.00013558*Z 

101 

SES JO=W/SQRT  <A)  *003 ( V) 

111 

RETURN 

112 

END 

6 

FUNCTION  FU(I,A,c) 
COMMON/ AUX/H,P,PK1,PK2, 

6 

X=A 

7 

Y=B 

10 

IF(A*9)5,10,5 

11 

10 

FU=0.0 

12 

RETURN 

13 

5 

SUM=SIMF(I,0. 0,5.0) 

20 

ER=0.01 

21 

DEL  =5. 0 

23 

20 

UP=OEL+5.0 

25 

A0DL=SIMF(I,DEL,UP) 

32 

DEL  =UP 

33 

TEST=ABS{AODL/SUM) 

36 

SUM=SUM+AODL 

37 

IF(TEST-ER) 15,20,20 

41 

15 

FU=SQRT (X*Y)*SUM 

47 

RETURN 

47 

ENO 
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SUBROUTINE  CONST (I) 

3  COMMON/  /SUX/H,R,PK1,PK2»BMU,X,  Y 

3  PR1=0.29 

5  PR2=0.2g 

6  PK1=SQRT  ((l.-2.*PRl)  /C2.^C1.-PR1)>  ) 

15  PK2=SQRT(C1.-2.»PR2) /(  c.*(l.-PR2))) 

24  READ  1,P 

31  1  FORMAT(F10.5) 

31  HH=0.1 

33  HH=10.0- 

34  HH=5.0 

36  HH=4.0 

37  HH=1.0 

41  HH=0.5 

42  HH=2.0 

44  H=1./HH 

45  PRINT  2,eMU,FRl,PR2,HH,P 

62  2  F0RMAT</////5X,*  MU2/MU1  =*F6.2,*  NUl  =*F4.2,*  NU2  =-*F4. 2///5X,*  A 

1/H  =^F4.2**  C21/PA  =»F4.2) 

62  RETURN 

63  END 


FUNCTION  Z(I,S) 

COMMON/ AUX/H,P,PK1, PK 2, €MU,X,Y 

PP=P*P 

C1  =  PK1-*PK1 

C2=PK2^PK2 

CC= 1 • “C  1 

GA=SQRT (S*S+C1/PP) 

G3=SQRT (S»S+1./PP) 

GC=SQRT {S*S+C2/8MU/PP) 

GD=SQRT (S*S+1./BMU/FF) 

AA=S*S+l./PP/2. 

AB=1.-BMU 
AC=S»S-GC*GD 
AD=fGB-G0)/AC/PP/2.*BMU 
AE= (G3+GD)/AC/PP/2.»BMU 
AF= (S^S-GA+GD) /AC/PP/2. ♦EMU 
AG=  (S*S+GA-‘G0)/AC/PP/2.  ♦EMU 
AH= (S*S-GB*GC)/AC/PF/2.*aMU 
AI= (S*S+Ga*GC)/AC/PP/2.*EMU 
AJ= (GA-GC) /AC/PP/2.^BMU 
AK= (GA+GC) /AC/PP/2.*BHU 
A1=-(A8*GB-A0) 

A2=AB»GE-AE 

A3=AA-BMU*S»S-AF 

A4=AA-BMU*S*S-AG 

A5=-AA+eMU*S»S+AH 

A6=-AA+3MU»S»S+AI 

A7=S^(AE^GA-AJ) 

A8=-S*< AE*GA-AK) 

BA=A1*A6-A2»A5 
BB=A3*A6-S’A2*A7 
EC=A4»A6-S»A2^A8  , 

BD=S^A1^A7-A3»A5 

a£=S*Al*A8-A4»A5 

ai=88/9A 

B2=BC/0A 

a3=BD/BA 

B4=e£/BA 

EA=2.*GA^H 

EB=2.*GE*H 

£C=(EA+£E>/2. 

£0=2, ^£0 

E1=EXP(-EA) 

E2=EXP(-EB) 

E3=EXP(-EC) 

E4=EXP(-£0) 

DL=B2+B3*E4+64*£2+ei^El 

D1=2.^PP/CC/GA/0L 

D2=AA^AA-S^S^GA^GB 

D3=B2-B3»£4 

D4=2.*AA^<Ge*(31^B4-B2*E2) -S*S»GA)*E3 
D5=  (AA^  AA'»-S*S^GA*GE)  ♦  (B4*£2-31^E1) 
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f-Ql* (02^03 ♦04+05) 

Z=CF-S) ♦eESJ0(S+X)*BESJ0(S+Y) 

RETURN 

END 


THIS^PROgIaM^EVALUATES’thFcoEFFICIENTS  FOR  SERIES 
OF  JACOBI  POLYNOMIALS  WHICH  REPRESENTS  A  LAPLACE 
INVERSION  -INTEGRAL  - 


INVERSION 

DIMENSION  A{5a),GLAM(50),PHI(50)  ,0(4,50) 
DIMENSION  BK  (101) ,TT (101) 
C0MM0N/2/TI,TF,0T,MN, eK,TT 
READ  l,NN,Mr4,MM 
F0RMAT(3I2)^ 

READ  2,TI,TF,0T 
FORMAT(3F10. 5) 

PRINT  99 

CALL^ SPL ICE ( GLAM , PHI , MM , C ) 

PRINT  101  -  - 

F0RMAT(/////5X,+  GLAM  PHI 

PRINT  102, (GLAH(I) ,PHI(I) ,1=1, MM) 
FORMAT(5X,F10.5,5X,F10.5) 

M11=MM-1  - 

FORMAT? ////5X,*  C(1,J)  C(2,J) 


A (50) ,GLAM(5Q) 


C(2,  J) 


C(3,  J) 


^^RINT  103, X (CC1,J) ,1=1,4) ,J=1, Mil) 

;  FORMAT(5X,F10.5,5X,F10.5,SX,F10 .5,5X,F10.5) 
PRINT  99 
00  10  I=1,NN 
READ  3, BET, DEL 
J  FORMATt2F10,5) 

I  FORMAT? />///§xf+aETA  =+F5.3,*  DELTA  =+F5.3) 
DO  11  L=1,MN 
AL  =  L 

S=1./(AL+3ET)/0£L  ,  ^  ^ 

CALL  SPLINE (GLAM, FHI, MM, C,S,G) 

F=G*S 

IF(AL-2.) 81, 82,83 
.  A(1)=(1.+0£T)+O£L*F 

I  A? 2)s(( l•♦SET)♦0£L♦F-A (1) )* (3.  + SET) 

GO  TO  11 
!  CONTINUE 
TOP=l. 

L1=L-1 

AL1=L1 

DO  12  J=1,L1 
AJ=J 

TOP=AJ=TCP 
»  CONTINUE 
L2=2»L-1 
eoT=i. 

DO  13  J=L,L2 
AJ=J 

aOT  =  (AJ+eET)  *BOT 
I  CONTINUE 
MUL=80T/T0P 
SUM^O.O 
DO  14  N=1,L1 
AN=N 

IF(AN-2.) 85,86,87 
?  T0D=1. 

GO  TO  88 
j  T0C=AL1 
60  TO  88 
r  CONTINUE 
TOD=l. 

ICH=L1- (N-2) 

DO  15  J=ICH,L1 
AJ  =  J 

TOO=AJ+TCO 
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15  CONTINUE 
88  CONTINUE 

eoc=i. 

JA=L1+N 
DO  16  J=L,JA 
AJ=  J 

eOO=BOD*(AJ+e£T) 

16  CONTINUE 
CO=TQD/EOD 
SUM=SUM+CO*A(N) 

14  CONTINUE 

ACL)=:MUL»(0£L*F-SUM) 

11  CONTINUE 

CALL  JACSERCD£L,A,B£T) 

CALL  NAMPLT 

CALL  QIKSET(e.O«0.O}O.Qt6.O)O.O»Q.O) 
CALL  QIKSAX(3,3) 

CALL  QIKPLT (TT»BK,101) 

CALL  ENDFLT 
10  CONTINUE 
999  CONTINUE 
RETURN 
END 
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6 
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SUBROUTINE  JACSER CD, C , 8) 

OIMENSICh  C(50)  ,SF  (50)  ,P(50) 

DIMENSION  BKClOl) ,TT(101) 

COMMON/2/TI ,TF,OT,MN,EK,TT 

TT(1)=0.0 

eK(l)=0.0 

LM  =  1 
T=TI 

12  T=T+DT 

X=  2,»EXP(-D*T) -1. 

CALL  JACOBI (MN,X,B,P) 

SF(1)=CC1)*P(1) 

DO  10  L=2,MN 

L1=L-1 

AL  =  L 

SF  (L)=SF(L1)  ♦C{L)*F{L) 

10  CONTINUE 
PRINT  97»T,X 

97  FORMAT( /////5X,*  T  =»F6.3,*  X  =^F10.5) 

PRINT  96 

96  FORMAT(///5X,*  I  C(I)  *,5X,*  N  F(T1  *) 
DO  11  1=1,6 

PRINT  95, I, C  (I)  ,I,SF(I) 

95  FORMAT(5X,I2,F10, 2,5X,I2,F10.5) 

11  CONTINUE 
LM=LM+1 
BK{LM)=SF(5) 

TT(LM)=T 

IFd.LE.TF)  GO  TO  12 

RETURN 

END 


SUBROUTINE  JACOBI  (N,  X,B  ,P8) 

C  THIS  PROGRAM  CALCULATES  JACOBI  POLYNOMIALS  OF  ORDER 

C  K-1  WITH  ARG  X  AND  PARAMETER  B  GT  -1 

7  DIMENSION  PB(N) 

7  AN  =  N 


10 

IF {AN-2.) 1,2,3 

12 

1 

PB(1)=1. 

14 

RETURN 

14 

2 

PB (1)=1. 

16 

pa(2)=x-e^<i.-x)  /2 

21 

RETURN 

22 

3 

aSQ=8+B 

23 

80NE=9+ 1. 

25 

PB(1)=1. 

26 

P8(2)=X-B*(l.-X)/2 

31 

DO  4  K=3,N 

33 

AK  =  K 
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AK1=AK-1. 

AK2=AK-2. 

K1=K-1 

K2=K-2 

CO  1=(  (2.*AK1) +B)  ♦X 

C01=( f2.^AK2) +B) *001 

G01  =  U2.^AK2)+B0N£)»(C01-eSQ)  ^ 

C02=2.*AK2»  < AK2+B)»( C2.*AK1) fB) 

C0=2.*AKl*(AKH-e)*((2.*AK2) +S) 

4  PB(K)=(C01*PB(K1> -C02*PB (K2) )/C0 

RETURN  - 

END 


SUBROUTINE  SPLINE  (X, Y , M »C ♦ XINT , YINT ) 

DIMENSION  X(50)  ,  Y(50)  »C(4,50) 

IFCXINT-Xd)  )1»10,11 

10  YINTsY(l) 

RETURN 

11  CONTINUE 

IFCX(M) -XINT) ltl2,13 

12  YINT=Y(M) 

RETURN 

13  CONTINUE 
K=M/2 
NsM 

2  CONTINUE 

IF(X<K) -XINT)3»14,5 

14  YINTsY(K) 

RETURN 

3  CONTINUE  ,  ,,  , 

IF(XINT-X  (K+D)  4»15t7 

15  YINT=YCK+1) 

RETURN 

YINTi(X^K+l)-XlNT)*{C(liK)*  (X(K+1)-XINT)**2>C(3,K>) 
YINT=YINT4-(XINT-X<K)  )*(d(2,K)*CXINT-X(K))»*2<-C(4,Kn 
RETURN 

5  CONTINUE  , 

IF  (X  (K-D-XINT)  b»  16»  17 

6  K=K-1 
GO  TO  4 

16  YINT=YCK-1) 

RETURN 

17  N=K 
K=K/2 
GO  TO  2 

7  LL  =  K 

K=  (N+K)  /2 
3  CONTINUE 

IF(X(K)-XINT)3, 14,18 

18  CONTINUE  ,  „ 

IF{XCK-l)-XINT)6,l6,ig 

19  N=K 
K=(LL+K)/2 
GO  TO  8 

1  PRI NT  101 

101  FORMATC*  OUT  OF  RANGE  FOR  INTERPOLATION  ♦) 

STOP 

END 


gSPo!!S!^oi:?f{i5?:5’(io).P<50,,E<50).C(4.50) 

DIMENSION  A(50,3)  ,6(50)  ,Z(50) 

MM=M-1 
GO  2  K=1,MM 
DCK)=X(K+1)  -X(K) 

P(K)=D(K)/6. 

2  E(K)  =  (Y<K  +  1) -YCK)  )/D<K> 

DO  3  K=2,MM 

3  e(K)=E(K)-E (K-1) 

A(l,2)=-l.-0(l)/0(2) 

A(l,3)  =  0  (1)  /D(2) 

A(2,3)=F(2>-P(1)*A(1,3) 
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A(2,2)  =  2.*(P(1)+F(2))-P(1)*A(1,2) 
A(2,3)=A(2,3)/A(2,2) 

8(2)=B(2)/A (2,2) 

DO  4  K=3,MM 

ACK,2)  =  2.»(P  CK-1)+P(K)) -F  CK-1)*A{K-1, 
B(KJ=B(K)-P(K-1)»B(K-1) 

- A(K,3)  =  P(K)  /A(K,  2) 

4  8(K) =a( K) /A (Kj2) 

Q=0  {M-2)/0  (M-1) 

A(M,1)=1.+Q+A(M-2,3) 

-  A(M,2)=-Q-A(M,1) ♦A(M-1,3) 

B{M)=BCM-2)-A(M,1)*B(«-1) 

Z(M)=B{M)/A(N,2) 

MN=H-2 


6 


7 


DO  6  1=1,  MN 

Z(K)=8(K)-A(K,3)*Z(K+1) 

Z( 1)=-A(1,2)»Z(2)-A(1,3)*Z(3) 
DO  7  K=1,MM 
Q=1./(6.^D (K)) 


C( 1,K)=Z(K)*Q 
C(2,K)=Z(K+1)*Q 


C(3,K)=Y(K) /0(K)-Z(K)^P<K) 
C(4,K)=Y(K*1)/0(K)-Z(K+1)*P(K) 


RETURN 


END 


3) 
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141 
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PROGRAM  3cTA { INPUT  , OUTPUT, PUNCH, PLOT »TAP£  99=PL0T) 
REAL  NON  (4)  ,F(4,4,1)  ,G(4,4)  ,0{*+)  ,PT(4) 

REAL  B{L ) ,C (4) 

REAL  LP(50) ,OTA(5Q) 
equivalence  (NON, 3) 

COMMON  K1,K2,<3,K4 
C OMMON/A UX/H ,P, PKl  ,PK2,BMU,X ,Y 
LP(1)=0, 0 
DTA(1)=G .0 
READ  2,K1,K2,K3,K4 
F0RMAT(I2) 

=  ORDER  OF  SYSTEM  OF  EQUATIONS 
=  NO.  OF  DISTINCT  KERNELS 
=  NO.  OF  DATA  POINTS  ^  , 

=  NO.  OF  DATA  SETS  TO  BE  EVALUATED 
UP  DATA  POINTS 
AK=K3 

DO  5  N=1,K3 
AN=N 

PT ( N)=AN/AK 
UP  INTEGRATION  MATRIX 
M=K3“2 
N=K3-1 
A=K3 

A=l./(3. ♦A) 

DO  1C  K=2,M,2 
D(K)=2.»A 
DO  15  K=1,N,2 
3  (0=4. »  A 
D  <K3)=A 

CULATE  NONHOMOGENEOUS  TERMS 

RHS=1.G 

DO  22  1  =  1, K2 

PRINT  9 

FORMATCIHI) 

READ  61,8MU 
FORMAT (FIG. 5) 

DO  999  11=1, K4 

00  35  N=1,K3 

NOra  N)  =RHS’SQRT  (PT(N)  ) 

CULATE  KERNEL  MATRICES 
CALL  CONST (I) 

DO  20  N=1,K3 
DO  2C  M=1,K3 
IF( M-N)25, 30 , 30 
F (M, N,I) =F ( N , M, I) 

GO  TO  2C 

F (M,N,I) =FU {I,PT( M),PT(N) ) 

CONTINUE 

CALL  CHANGE {F,G,0, I) 

CAlL  LI  N  EQ  (  G  ,  i ,  C  ,  Kj) 

DO  40  L=i,K3 
PRINT  6, PT (L) ,NON (L) 

F0RMAT{5X,F8.4,  Fi5.6) 

C  ONTINUE 

LP(II^l) =N0N(<3) 

DTA (II  +  l  )  =P 
CONTINUE 

PUNCH  66  ,  (OTA  (IX)  ,LP(IX)  ,IX=1,  19  ) 

FORMAT(2F10 .5) 

CALL  LAPINV(OTA,LP) 

CONTINUE 

END 


FUNCTION  SIMPCI,A,8) 
COMMON/AUX/H.P,  PK 1  ,PK  2  ,  B  MU  ,  X  ,  Y 
DEL=0.25* (3-A) 

I F(DEL)LG,45, 50 
45  SIf1P  =  0.C 
RETURN 
50  CONTINUE 

3A=Z(I,A  )4-Z  (1,3) 
SB=Z(I,A+2.»0EL) 
3C=Z{I,AtCEL)i-Z(I,A  +  3.»0£L)  • 
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51=(DEL/3.)*{5fl+2.’‘S3+4.»3C) 

IF(S1.EQ.0.  Q),  50  TO  45 
K  =d 

35  33=SB+SC 

0EL=0.5»OEL 
SC=Z  (I,A«-0£LJ 
J  =  K-1 

00  5  N=3,J,2 
AN=N 

5  3  C=SC«-Z(I,Ai-AN*DEL) 

3 2= (DEL/ 3.)  ♦(SA+2.»Sa+4.»SC) 

0IF=A3S( (S2-S1) /31) 

ER=C.C1 

IF(OIF-£R) 30, 25,25 
30  SIMP=S2 
RETURN 
25  K=2»K 

31=3  2 

IF  (K-204d  )  35, 35,  0 

40  PRINT  42,1, A, 3 

42  F0RMAT(5X,»  INT.  DOES  NOT  CONVERGE  ^,I3,2F9.4) 
PRINT  60  ,X,  Y 
60  FORMAT(2Fia .5) 

DO  70  J=l,iO 
0IP=U 

OIF=DIP/10. 

.'I  =  Z(I,DIP) 

PRINT  6C,W 
70  CONTINUE 
call  EXIT 
END 


SUBROUTINE  CHANGE ( F,G,0, I) 
7  REAL  F(4,4,  1)  ,G (4 ,4) ,0(4) 

7  COMMON  K1,K2,K3,K4 

7  DO  10  N=1,K3 

10  DO  1C  M=1,K3 

11  G (M, N)  =F (N,N, I) »D(N) 

24  1C  CONTINUE 

3C  CO  20  N=1,K3 

31  2C  5  (N,N)=G  (N,N)  <-1.3 
4C  RETURN 

41  E  iO 
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SUBROUTINE  L I NE Q ( A , 3 , T, N ) 
REAL  A (N, N) , 3 (N ) , T (N) 

DC  5  1=2, N 

5  A  (I,1)=A  (I,l)/A(i,l) 

DO  10  K=2,N 
M=K-1 

00  15  1  =  1, N 
15  T(I)=A(I,K) 

DO  20  J=1,M 
A  (J,  K)=T  (J) 

J1=J+1 

DO  2C  I=Jl,N 
T  (I)=T(I)-A(I,J)»A{J,!<) 

20  CONTINUE 

A (K, K)=T (K) 

IF(K.EQ.N)  GO  TO  10 
M=K<-1 

DO  25  I=M,N 

25  A (I, K)=T (I) /A (<, K) 

10  CONTINUE 
*  3AGK  SUBSTITUTE 
DO  30  1  =  1, N 
T  (I)=B(I  ) 

N=H-1 

IF(N,GT. N)  GO  TO  30 
DO  30  J=M,N 
B (J) =3{ J) -A { J,I ) ♦T  (I) 

30  CONTINUE 

DO  35  1  =  1, N 
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137  K=N+1-I 

141  3 (K>=T{<)/A (K,<) 

146  K1=K-1 

15a  IF(K1.£Q.4)  GO  TO  35 

151  DO- 35  J1=1,K1 

152  J=K-J1 

154  T {J)=T(J)-A{J»<)»aCK) 

162  35  CONTINUE 

167  RETURN 

16  7  £  NO 


3 

5 

7 

12 

13 

15 

16 
20 
22 

24 

25 
27 
31 
34 
34 
36 
41 

44 

45 
47 

52 

53 
55 
57 
61 
63 
65 
6  7 
71 
73 
75 
77 

IDl 

111 

112 


INUNCTION  BESJO(A) 

IF(A-3.} 5,5,10 
5  3  =A*A/9. 

W=1.-2.2499997»B 
Z  =  B»B 

W=W+1.2656203*Z 

Z=Z*8 

W=W-.3163866*Z 
Z  =  Z*B 

W=W+ .0444479*2 

z=z*a 

W=W-.0C39444*Z 
Z  =Z*B 

BESJC=W+. JG021*Z 
RETURN 
IQ  B  =  3./A 

W=. 79786456-. 00000077*8 
'7=A-.78539816-.  04166397*3 
Z=3*3 

W=W-, 0055274*2 
V=V-.0CGC3954*Z 

z=z*a 

W=W-.00G09512*Z 
7=^+.00262573*Z 
Z  =Z*  B 

W=W+.QC137237*Z 
\/  =  \/-.aOO  54i25*z 
Z  =  Z*B 

W  =N-  .  uC  C  72  3  0  5  *Z 
y=7-.QGC29333*Z 
Z  =Z*B 

N=W+,0tiCl4476*Z 
\/=V<-.00C  13558*Z 
3ESJ0=W/SQRT  ( A)  *COS(\/) 
RETURN 
E  NO 


FUNCTION 


6 

OONMON/AUX/H,P,PK1 , 

6 

X=A 

7 

Y=8 

1C 

IF(A*a)5,lC ,5 

11 

10 

FU=0 .0 

12 

RETURN 

13 

5 

3UM=SIMP (I,  0 . 0, 5.  a ) 

20 

£R=a .01 

21 

DEL  =5.C 

23 

20 

Uo=OEL+5 . 0 

25 

AODL=SIMP(I,OEL,UP) 

32 

OEL  =UP 

33 

T£ST=A3S (AOOL/SUM ) 

36 

3U.'1=SUH4  AOOL 

37 

IF  (TEST-  ER)  15,20, 20 

41 

15 

FU=SQRT (X*Y ) *3UM 

47 

RETURN 

47 

E  NO 
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3U3R0UTI N£  CONST ( I  ) 

3  C0NM0N/AUX/H,P,P<l,PK2,8hU,X,Y 

3  PR1=0.29 

5  PR2=C.29 

6  PK1  =  SQRT  ((l.-2.*PRl)/{2.*{l.-PRin) 

15  PK2=SQRT ((i.-2.*PR2)/(2.*(l.-PR2))) 

24  READ  1,P 

21  1  FORMAKFIO.  5) 

31  HH=0.1 

33  HH=1C.C 

34  HF=5.0 

36  HH=4.0 

37  HH=0.5 

41  HH=1.G 

42  HH=2.0 

44  H=1./HH 

45  PRINT  2»BMU ,PRl,PR2tHH,P 

62  2  FORMAT(/////5X»  »  MU2/HU1  =^F6.2,*  NU  1  =^F4.2,*  NU2  =  *F4 . 2/ // 5X ,  ♦ 

1/H  =*F4.2,*  C21/PA  =»F4.2) 

62  RETURN 

63  E  NO 
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FUNCTION  Z(I,S) 

COMMON/AUX/H  ,P,  PK 1  »PK 2,  BMU  ,  X  ,  Y 

PP=P*P 

Zi=PKl*FKl 

C2=PK2*PK2 

SA=SQRT(S»5+Ci/PP) 

3B=SQRT{S»S+1./PP) 

GC=SQRT(S*3+C2/3MU/PP) 

GO=SQRT(S*S+l ./BMU/PP) 

AA  =  3*S+1  ./PP/2. 

AB=i  ,-BKU 
AC=3»S-GC^GO 
AD=(G8-GD)/AC/PP/2.»flMU 
AE= (GB+GO) / AC/PP/2 .♦BMU 
uF=  (S'^S-GA^GD)  7AC/PP/2.*B:1'J 
AG=  (S*Si-GA^GD)/AC/FF/2.»a,4U 
AH=  (S'^S-GB'^GO/AC/PP/Z.^BMU 
AI= (S*S+G3»GC)/AC/PP/2.*BMU 
A j= (GA-&C)/AC/PP/2.*8MU 
AK= (GAtGC)/AC/PP/2.»3RU 
A  1  =  -  (A3'-G3-AD) 

A2=AB*GB-AE 
A3=AA-BNU*S»S-AF 
A4=AA-BNU*S*S-AG 
A5  =  -AA«-BMU»S^S^AH 
A6=- AA+BMU^S»5+AI 
A7=S»(AB^GA-AJ) 

A8=-5*(A3*GA-AK) 

3A=A1*A6-A2»A5 
3a=A3»A6-5^A2*A7 
3G=A4^Ab-S^ A2*AS 
3D  =  S^A1*  A7- A3*A5 
3E=S ♦A1*A3-A4*A5 
3i=3B/EA 
3  2  =  3  C/ BA 
3  5=3  0/3  A 
34=BE/BA 
EA=2.*GA*H 
£B=2.*&B*H 
E  C=(EAi-EB)/2. 

ED  =  2  .*EC 
E1  =  £XP(-EA) 

E  2  =  EXP (“EB) 

E3=EXP(-£C) 

E4=EXP  (-EO) 

DL  =  B2  +  33'^E^-BH»E2-ai^El 
Dl=2  .♦PP/CC/G3/DL 
C)2=AA’^AA-S’^S»GA*G3 
33  =  32-33»E-+ 

D4=2.'^AA*(G5*{3i»B4-B2*B3)  -S’^S’^GA)  ^£3 
05=:(AA-^AA  +  S*S*GA^‘GB)*(B4^E2-B1*£1) 
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RETURN 

END 
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INVERSION, INTEGRAL 

OIMENSICN  A15?'jf'-4H|501j'>Hll5(l>.C(4.5t)) 

DIMENSION  3K(101) .TTtlCl) 

C0MM0N/2/TI»Tr tOT ,KN»BK,TT 
READ  1,NN,MN,MM 

1  F0RMATt3I2l 
READ  2,TI,TF»DT 

2  FORMATCSFIO .5 ) 

PRINT  99 

101  FORmIt(/////5X,  ♦  ohr^/Ti  T-i  MM) 

PRINT  10  2  f  (GLAM  ( I  )  t  PHI  ( I )  *  I- 1»  NM  ) 

10 2  FORM AT (5 X f FIO .5 » 5 X »F10. 5 ) 

M11=MM-1 

300  co^maT(////5X»*  C(1*J)  C{2tJ)  ^(3 

3RINT  99 
DO  10  I=1»NN 
REAO  3, bet, del 

3  FORHAT(2Flg.5) 

99  FORMATC///7/5x7*3£TA  =*F5.3,*  DELTA  =*F5.3) 

DO  11  L=1»MN 
A  L= L 

C  ALu'^SPLINllGLAMt  PHI  ,  MM,  C,  S  »  G) 

F=G*S 

IF(AL-2.) 81,02,83 
di  A  {1)  =  (1.  ♦•BcT)  ♦Dc.L*F 

32  A  (2l=(i2  ,+aET)  ♦OEL»F-A(l)  )  *  (  3.  *-BET) 

GO  TO  11 

33  CONTINUE 
TOP=l. 

Li=L-l 

ALl^Ll 

DO  12  J=1»L1 

AJ=J 

top=aj»top 

12  CONTINUE 
L2=2*L-1 
BOT=l. 

3  C  13  J  =  L  »L2 

30T=  ( AJ«-3ET  )  *33T 

13  CONTINUE 
MUL=BOT/TOP 
3UM=C.C 

00  14  N=1,L1 

an=n  ^ 

I F( AN-2.  )  85,80,87 
85  T00=1.  , 

SO  TO'88 
do  T0D=AL1 
GO  TO  88 
37  CONTINUE 
TOD=l. 

ICH=L1- (N-2) 

•DO  15  J=ICH,L1 
•AJ=J 

T0:=AJ*T00 
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250 

15 

C  ONT INUE 

252 

83 

CONTINUE 

252 

3  0  0=  1 . 

254 

J  A=i_l*-N 

256 

DO  16  J=L,JA 

260 

A  J=  J 

261 

300=B0D»  (Aj^-BET  ) 

254 

16 

CONTINUE 

266 

CO=TOO/BOO 

270 

SUM=SUH+CO»A (N) 

273 

14 

CONTINUE 

275 

A (L)=MUL»{OEL*F-SUM) 

301 

11 

CO(4TINUE 

334 

CALL  JACSER (DEL, A ,8ET) 

306 

CALL  NAMPlT 

307 

CALL  QIKSET (6. C, 0.0, 0.0, 6. 0,0. 0,0.0) 

3  13 

CALL  QIXSAX(3,3) 

315 

CALL  QIKPLT (Tr,BK,iai) 

320 

call  endplt 

321 

IQ 

CONTINUE 

325 

999 

CONTINUE 

325 

RETURN 

3  26 

END 

6 

6 

fc 

6 

7 

IG 

11 

12 

14 

24 

26 

32 

33 
35 
3o 
43 


45 

55 
55 
6 1 
61 


65 
135 
105 
111 
113 
115 
117 
121 
1  22 


SUBROUTINE  J ACS ER ( Ot C ,B ) 

DIMENSION  C(5G) ,SF(5Q)»P{50) 

DIMENSION  BK{101} ,11(101) 

SOHHON/2/TI  ,TF,  CT  ,NN,3K,TT 

TT(l)=C.O 

3K(1) =2. 0 

LM=1 

T=TI 

12  T=T+DT 

X=2.*EXF (-0*T)-1. 

CAuL  JACOSI (HN, X, 3,P) 

SF(1)=C{1>*P(1) 

DO  IG  L=2,MN 
L1=L-1 
A  L  =  L 

SF  (L  )  =SF  (LI  )  (  L  )  ♦P  (L  ) 

10  SOIJTINUE 
=RINT  97,T,X 

97  FORMAT(/////5X, ♦  T  =»F6.3,*  X  =»F12.5) 

=RINT  96 

96  FORMAT(///5X,»  I  Cd)  *,5X,*  N  F(T)  ^) 
uO  11  1=1,6 

ORINT  95  ,I,SF{I) 

95  FORMAT(5X,I2,F10.2,5X,I2,F1C.5) 

11  CONTINUE 
LM=Lh+l 
3K{LM)=SF (5 ) 

TT(lN)=T 

IF(T.LE.TF)  GO  TO  12 

RETURN 

END 


7 

7 


SUBROUTINE  J A CO BI { N, X , 3 , P3 ) 

C  THIS  PROGRAM  CALCULATES  JACOBI  POLYNOMIALS  OF 

C  <-l  WITH  ARG  X  AND  PARAMETER  B  GT  -1 

DIMENSION  P3(N) 

A  N=  N 


iC 

IF (AN-2. ) 1,2,3 

12 

1 

PB (1)=1. 

14 

RETURN 

:4 

2 

P3(l)=l. 

PB{2)=X-8^(l.-X)/2 

RETURN 

;2 

3 

3SO=B^3 

13 

30NE=B+1 . 

cO 

PB(1)=1. 

2b 

=  3(2)=X-3^‘(l.-X)/2 

::1 

30  4  K=3,N 

23 

A  K=K 

ORDER 
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34 

36 

40 

42 

43 
46 
51 
56 
64 
71 

102 

1Q3 


aKl= AK-1. 
a<2=AK-2 . 

K1=K-1 

— 2 

301=  (  C2  .  *AK1)  ^5  ) 

301=  !  (it  *AKi)ipNE^^(C01-aSQ) 

4  P3{K)=(C01»PB(K1) -C02*PB(<2» )/C0 
RETURN 
■  ENO 


11 

11 

13 

14 

15 
15 
21 
23 
23 
23 

25 

26 
26 
32 

34 

35 
35 
41 
43 
43 
43 
54 
65 
65 
65 
7C 
72 
72 

74 

75 
77 

100 
IOC 
102 
103 
103 
106 
10  6 
111 

113 

114 

115 
121 
121 
123 


S  U3ROUTI NE  SP4lNECX»YfM*C»XINT»YINT) 
DIMENSION  X  (  5  0)  ♦  Y  (  53)  » C  ( 4»  50  ) 

IFIXINT-Xd)  )  !♦  10  »11 
13  YINT=Y(1) 

RETURN 

11  C0^4TINUE 

IF(X («)-XINT) l»12f 13 

12  YINT=Y(M) 

RETURN 

13  30NTINUE 
<=rV2 
N=M 

2  CONTINUE  ,  ^  _ 

IF(X  (K)-XINTI3»1‘»*5 

14  YINT=Y(K) 

RETURN  ^  , 

3  CONTINUE  ,  __ 

IF{  XINT-X(K<-1) )  4f  15»7 

15  YINT=Y(K  +  1) 

RETURN 

RETURN 

5  CONTINUE  . -  , 

IF(X  (<-l )-XlNT) 6» lc»l7 

6  K=K-1 
30  TO  4 

16  YINT=Y(K-1) 

RETURN 

17  N=< 

K=K/2 
GO  TO  2 

7  LL=K 

<= (N+K)/c 

<J  C  CNTINUE  a 

IF(X  (K)-XINT)  3t  18 

13  CONTINUE  ,,  ,Q 

IF(X  (K-D-XINT)  6  ♦  I6tl9 

19  N=< 

K=(LL+K)/2 
GO  TO  8 

loi  FORMsn*  OUT  OF  '?»NGE  FOR  INIERFOLAI  ION  M 

STOP 

ENO 


7 

7 

7 

11 

12 

15 

20 

26 

27 

34 

37 

Ll.  1 


§f.iSfl'oK%f55?^f!i5|:25io.5P<5D).E(50., 0(4,5=. 
DIMENSION  A(53»3)»8(50>*Z(50) 

MM=M-1 
DO  2  K=t»MM 
□  IK) =X(K  +  1) -X (K ) 

=  (K)='{Y(Kd)-Y(K)  )/0(K) 

iw=r7i5-£(K.i) 

a(l,2)=-l.-0(l»/D<2J 

A (lt3)=0 (1) /0(2 

A(2,3)=P(2)-P(i)*A(lt3) 


K)  ) 
*<)  ) 
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44 

5Q 

51 

53 

54 
61 
65 
7D 
74 
76 

101 
10  5 
112 
114 
116 
117 
1  20 
1  27 
133 
135 
140 
143 
146 
154 

4  rr 


a(2,2)  =  2.*(P(l)+-P{2))-Pll)''Atl»2) 

4  (2, 3)=A  (2, 3)  /4  (2 , 2) 

B (2)=B(2)/A(2,2) 

DO  4  K=5,iM 

4{K»2)=2.*(P(K-1)4P(K))-P(K-1)»A(K-J 
3  (K)=B(K)-P(K-1)»9(K-1) 

A (K,  3)=P (K) /A (<  ,2  ) 

3 (K)=3{K)/A(K,2) 

3=0 ( M-2) /O (M-i) 

A  (M,  l)=l.•^Q  +  AlM-2.3) 

A  (M, 2)  =  -Q-A (M,1)*A (M-l»3) 
3(M)=3(H-2)-A(H,l)*3{M-l) 

Z  {M)=a{M)/A(M,2) 

MN=M-2 
30  6  1=1, MN 
<  =  M-I 

Z{K)=3(K)-A(K,3)»Z  (K<-1) 
Z{1)=-A(1,2)*Z(2)-A(1,3)*Z(3) 

30  7  K=1,MM 
Q=l./ (6. *0(KI ) 

3(1,  K)=Z  (K)  *0 

C  (2,K)=Z  (K^-l)^a 

3  (3,K)  =  Y(  <)  /D(K)  -  Z  (K)  ♦P(K) 

C  14,  K)=Y  (K^H)  /O  (K)  -Z  (Ki-l  )»P  (K) 


RETURN 
£  NO 


,3) 
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